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ABSTRACT 

We present observations of six transits and six eclipses of the transiting planet system HD 189733 
taken with the Spitzer Space Telescope IRAC camera, at 8 microns, as well as a re-analysis of previously 
published data. We use several novel techniques in our data analysis, the most important of which is 
a new correction for the detector "ramp" variation with a double-exponential function which performs 
better and is a better physical model for this detector variation. Our main scientific findings are: (1) 
an upper limit on the variability of the day-side planet flux of 2.7% (68% confidence); (2) the most 
precise set of transit times measured for a transiting planet, with an average accuracy of 3 seconds; (3) 
a lack of transit-timing variations, excluding the presence of second planets in this system above 20% 
of the mass of Mars in low-order mean-motion resonance at 95% confidence; (4) a confirmation of the 
planet's phase variation, finding the night side is 64% as bright as the day side, as well as an upper 
limit on the night-side variability of 17% (68% confidence); (5) a better correction for stellar variability 
at 8 micron causing the phase function to peak 3.5 hours before secondary eclipse, confirming that 
the advection and radiation timescales are comparable at the 8 micron photosphere; (6) variation in 
the depth of transit, which possibly implies variations in the surface brightness of the portion of the 
star occulted by the planet, posing a fundamental limit on non-simultaneous multi- wavelength transit 
absorption measurements of planet atmospheres; (7) a measurement of the infrared limb-darkening 
of the star, which is in good agreement with stellar atmosphere models; (8) an offset in the times of 
secondary eclipse of 69 seconds, which is mostly accounted for by a 31 second light travel time delay 
and 33 second delay due to the shift of ingress and egress by the planet hot spot; this confirms that the 
phase variation is due to an offset hot spot on the planet; (9) a retraction of the claimed eccentricity 
of this system due to the offset of secondary eclipse, which is now just an upper limit; and (10) high 
precision measurements of the parameters of this system. These results were enabled by the exquisite 
photometric precision of the Spitzer IRAC camera; for repeat observations the scatter is less than 0.35 
mmag over the 590 day time scale of our observations after decorrelating with detector parameters. 
Subject headings: stars: planetary systems 



1. INTRODUCTION 



The planet system HD 189733 (jBouchv et al.l 120051 ) is 
one of the best studied transiting planet systems due to 
two factors: its close proximity to our Solar System, mak- 
ing its star one of the brightest transit host stars, and the 
large size of the planet relative to the star, making the 
transits particularly deep. After the secondary eclips e 
was first detected for this planet bv lDeming et al.l (J2006D . 
iKnutson et al.l ()2007| ) made a precise measurement of the 
phase variation of the planet over slightly more than half 
of an orbital pe riod using the 8 m icron Infrared Array 
Camera (IRAC, iFazio et all I2004D on the Spitzer Space 
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Telescope (jWerner et al.ll200"^ . I n addition to yielding a 
longitudinal map of the planet (iCowan and Agoill2008D 
which indicated an offset peak in brightness, attributed 
to a dvection of energy by a super-rotating equatorial 
iet (iShowman and GuillotI 120021 : ICooper and ShowmanI 
|2005[) ■ this observation also yielded the most precise mea- 
surement of the depth of secondary eclipse, as well the 
most precise times of transit and secondary eclipse, for 
any extrasolar planet. This motivated us to propose 
additional observations of six transits and six eclipses 
of this system with t he goals of looking fo r secondary 
eclipse variability (e.g. iRauscher et aI]|2007D . looking for 
trans it timing variations due to other planets in th e sys- 
tem (jAgol et all 120051 iHolman and Murray! I2005D . im- 
proving t he measurement o f the atmospheric absorp- 
tion (e.g. iTinetti et al.l 120071 1 . and improving the mea- 
sured system parameters for better charac terization of 
the planet, host-star, a i id orbit properti es (jWinn et al.l 
[20071: llbTres et al.l[200l IPont et al.l[2007l) . 

The favorable properties of HD 189733 have al- 
lowed detections of planet absorption and emis- 
sion features, yielding possible evidence for wa- 
ter, sodium, methane, carbon dioxide, and car- 
bon monoxide, as well as Rayleigh scattering at 
blue wavelengths (Grillmair ct al. 2007; ITinetti et al 



2007"; 'Barnes'et al."'2007; Rcdficld ct al. '2008; 'Barman 
2008; ,Swain et al.. .2008, ,2009; .Charbonneau ct a" 
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2008t iSing elaLl [200I iMadhusudhan and Seageil[200l 



Pont et al.ll2008l : iLecavelier Pes Etangs et al.ll200S 

Despite being an a ctive star (" Moutou et al.l 120071 : 
iHenrv and Winnll2008t ) which affects radial velocity mea- 
surements, the planet mass is measured precisely to be 
Mp = 1.13±0.03Mjupitcr, while the radial velocity mea- 
surements constra in the eccentricity to be e < 0.008 
(jBoisse et al.ll2009() . T he planet radius is s lightly larger 
than that of Jupiter jBakos et"al] l2006at iBaines et alj 
120071: iWinn et all 120071 ). The orbit of the planet is well 
aligned with the sp in axis of the star (jWinn et al.ll2006l : 
iTriaudet al1[2009[ i. 

The deluge of observational constraints on this system 
has inspired a wide range of theoretical modeling. In 
particular, the measured phase variation can be qual- 
ita tively explained b y gene ral cir culation models, such 
as iShowman et aP ()2009al ) and iRauscher and Menoul 
()2010t ). while two- hemisphere models have more diffi- 
culty exp laining the shape a nd variation of the phase 
function (iBurro ws et al.l l2008D : see the recent review by 
IShowman et al . (2009b) for a detailed discussion of these 
models. 

After describing our observations (ij2]), we give an ac- 
count of our preliminary data reduction (Sj3]), describing 
our outlier rejection f iJ3.1l) and choice of centroiding al- 
gorithm fi i3.2p . In section m we discuss aperture photom- 
etry, background subtraction (i j4.ip . and then detail our 
correction for the detector ramp variation f §4.2|) . includ- 
ing a new double-exponential model for the ramp (i i4.2.3p 
and its performance (i 34.2.4|) . We complete the descrip- 
tion of data reduction with a discussion of our choice 
of aperture size ( iJ4.3p . conversion to Barycentric Julian 
Date ( ij4.4p . and error analysis (S3]). 

With the preliminary fit for the ramp, we then simulta- 
neously fit to the stellar variability and planet variability, 
outside of eclipse or transit, and demonstrate the high 
precision of Spitzer IRAC (©. We then fit the photom- 
etry with transit and eclipse models {^, and show that 
the secondary eclipse depth offset can be explained by 
light-travel time and the offset hotspot ( §6.5|) . We com- 
pute a new ephemeris from our data ( §6.3|) . and use the 
times of transits to place limits on the presence of com- 
panion planets ( H6.4p . We show that the transit depth 
appears to vary, which we hypothesize is due to varia- 
tions in the stellar surface brightness within the path of 
the planet ( §6.6p and we show that the day-side planet 
flux measured from the secondary eclipses appears not 
to vary within the uncertainties ( §6.7p . We discuss these 
results and compare to models in the conclusions ([J7]). 

A prelimina r y ana lysis of these data were presented 
in lAgol et all (|2009f) : however, we have since made 
significant improvements in the analysis, in particular 
an improved ramp function, so the results presented 
here are more re liable. These data were also used by 
iCarter and WinnI (2010) to place a constraint on the 
oblateness of HD 189733b, while for the purposes of this 
paper we assume the planet to be spherical. 

2. OBSERVATIONS 

We were awarded Spitzer Guest Observing time during 
Cycle 4 to observe six transits and six secondary eclipses 
of HD 189733 with IRAC Channel 4 (PI: E. Agol, pro- 
gram ID 40238). For each visit we obtained 44,160 ex- 
posures of 0.4 second each over 5 hours each. We also 



re-ana lyze the transit and eclipse from iKnutson et al.l 
(|2007D for a total of seven eclipses and seven tran- 
sits. We chose IRAC channel 4 (8 ^m) as it has been 
demonstrated to b e the most stable IRAC band (e.g. 
ICowan et al.l 120071 and references therein). Due to the 
brightness of the host star we made the observations in 
sub-array mode; this mode allows shorter exposure times 
(0.4 sec) and faster readout, but sacrifices the larger field 
of view (32x32 pixels rather than 256x256 pixels, where 
each pixel is 1'.'2). We turned off dithering which is re- 
quired for high precision photometry due to the array- 
dependent sensitivity and detector ramp. We carried out 
aperture photometry with a range of radii from 1 to 7 
pixels in 1/2 pixel increments. 

3. DATA REDUCTION 

We performed our data reductions starting with the 
Basic Calibrated Data (BCD) processed with version 
S16.1.0 of the Spitzer /i?^C pipeline. These data are cor- 
rected for dark current, fiat field variations, and detector 
non-linearity; they are also converted to units of fiux in 
mega-Jansky per steradian (MJy/sr). After downloading 
and organizing the data, we first converted the images to 
units of photon counts (i.e. electrons) by multiplying by 
the gain (fits header keyword GAIN = 3.8 e~/DN) and 
exposure time (EXPTIME = 0.32 sec), and dividing 
by the flux conversion factor (FLUXCONV = 0.2021 
MJy/sr per DN/sec). For our 0.32 sec exposure time 
this amounts to multiplying each pixel by 6.01682 e~ per 
MJy/sr. The elapsed time per exposure is FRAMTIME 
= 0.4 seconds due to a 0.08 second readout. We did not 
apply corrections for variation in pixel area or corrections 
to the flat field for a stellar source; however, we did es- 
timate the impact of these corrections, and found them 
to be negligible. 

3.1. Outlier rejection 

After conversion to counts, we fiagged and cleaned the 
images of outliers, such as cosmic rays. We cleaned the 
images at the pixel level by rejecting outliers in the time 
series for each pixel. This worked well since the telescope 
pointed at nearly the same location for the entirety of 
each of our observations, so the flux in each pixel stayed 
relatively constant making outliers easy to flag. The pixel 
flux is affected by pointing variations, discussed in sec- 
tion 13.21 so we did the outlier rejection by taking the 
difference between the pixel time series and a 5 exposure 
(2-second) running median of the time series. The du- 
ration of the median was chosen to be shorter than the 
shortest of the pointing excursions, one of which occurred 
in the middle of the transit of the phase-function obser- 
vation with a 1-pixel pointing change over 4 seconds. 

For each pixel, the median-subtracted time series was 
sorted. Then, the standard deviation was computed from 
the 68.3% confidence limits on the median-subtracted 
time series. Finally, outliers were flagged in the median- 
subtracted time series that differed by more than A-a 
from zero. The flagged pixels were then replaced by the 
5-exposure median. 

This procedure performed well in eliminating cosmic 
rays and rogue pixels. After it was carried out, the pho- 
tometry showed no significant outliers and, as we show 
below, was very close to the photon-noise limit. 



Transits and Secondary Eclipses of HD 189733b 



3.2. Centroiding 

Spitzer is affected by pointing variations that cause the 
star to change position shghtly on the detector; this re- 
quires accurate centroiding to perform precise photome- 
try. There are four varieties of pointing variations we ob- 
serve in the data: small amplitude, short-timescale "jit- 
ter" which appears by eye as a damped random-walk; 
a periodic pointing fluctuation which occurs on a time 
scale of '^ 1 hour with an amplitude of about 0.1 pixel; 
gradual drifts over long timescales; and occasional short 
time scale sharp excursions, as mentioned in section [3. II 

In the inital stages of our data reduction, we found that 
the photometry was extremely sensitive to the pointing 
drifts. Since the 8 micron camera is undersampled, vari- 
ations in pointing at the ~ 0.1 pixel level can lead to ~ 
10% variations in the flux of the central pixels. For large 
photometric apertures, this is not a problem since small 
changes in the centroid do not change the enclosed flux 
much. However, large apertures contain pixels with lower 
illumination which have a longer-lasting detector ramp 
(see section 221) , so a smaller aperture containing higher 
illumination pixels is more desirable since these pixels 
have a ramp that saturates more quickly. Thus, we real- 
ized that a very accurate centroiding algorithm would be 
necessary for using smaller apertures, so we set out to test 
a wide range of different centroiding algorithms to see 
which performed best. We tried several centroiding al- 
gorit hms: flux- weighted cen troiding (e.g. IKnutson et al.l 
I2007f) . parabolic fitting (e.g. iTodorov et al.ll2010D. Point 
Response Function (PRF) fitting fe.g. 'Lau ghlin et ahl 
|2009), and 2-D Gaussian- fitting (e.g. Desert e t al.ll2009[ ). 
The simplest algorithms to implement are the first two 
since these involve no optimization; the last two involve 
iterative non-linear optimization, but the PRF fitting 
turned out to be too slow and problematic to implement. 

We first tested the centroiding algorithms by creat- 
ing simulated jitter using the IRAC channel 4 PRFs. 
From these tests we found that the 2-D Gaussian fitting 
gave the least scatter in the derived centroid relative to 
the input centroid; we found that keeping the x and y 
standard deviations the same gave as good a fit to the 
centroid as allowing the two to vary independently. We 
next ran the algorith ms on the two stars ( target star and 
M-dwarf companion; 'Bakos et al.' 2006b) in the phase- 
function data from Knutson et al. (2007) and on the two 
stars in the observations of HD 80606 (jLaughlin et al.l 
120091 ). These tests were critical since these were long 
time series so the stars had time to drift a significant 
fraction of a pixel, and the data contain noise. We first 
computed the centroid of both stars in each image {xi,yi 
and X2,y2), then subtracted the x and y coordinates for 
each pair of stars (a;i — a;2 and j/i — 2/2), and finally binned 
the x and y differences until the standard deviation of 
the X and y differences reached a minimum. A perfect 
centroiding algorithm ought to have perfect tracking be- 
tween the two stars, resulting in a standard deviation 
due only to the photon shot noise and finite spatial res- 
olution of the instrument. Various choices can be made 
for each of these algorithms, such as what portion of the 
array to fit or whether to smooth the data first before 
centroiding, so we spent some time experimenting with 
these and other choices. 

In short, we found that the 2-D Gaussian performed 



the best of all the centroiding algorithms. The algorithm 
selects a 7 X 7 sub-array from the image centered on the 
brightest pixel of a star. It then fits a 2-D Gaussian to 
this sub-array, allowing the center (centroid), amplitude, 
and width to vary; four free parameters in all for each 
image. We used the mpcurvef it .pro routine which im- 
plements a non-linear Leven berg-Marquardt al gorithm to 
optimize these parameters ()Markwardtl 120091) . For the 
HD 189733 phase-function observations, the scatter in 
the data were 0.0018 pixels in x and 0.0051 pixels in y 
when binned by 512 exposures (205 seconds), while for 
HD 80606 the scatter was 0.0015 and 0.0021 in x and y 
when binned by 4 exposures (56 seconds); further bin- 
ning resulted in minimal decrease of the scatter. The 
second best centroiding algorithm was the flux-weighted 
algorithm which had standard deviations of 0.016 and 
0.032 in X and y for HD 189733, and 0.0019 and 0.011 
in X and y for HD 80606. Thus, the 2-D Gaussian cen- 
troid performed better by a factor of ^ 5 than the flux- 
weighted centroid; not only that, but the scatter in the 
flux-weighted centroid is due to a systematic error, while 
the scatter in the 2-D Gaussian centroid is almost com- 
pletely random. This can be seen in Figure[T]which shows 
that the 2-D Gaussian centroid has weaker correlation of 
centroid difference of the two stars versus the centroid of 
one of the stars; on the other hand the other centroid- 
ing techniques show signficant correlation between the 
offsets of the two stars and the pixel position of one of 
the stars, indicating a systematic error in the centroid 
determination. 

Applying these centroiding algorithms to our twelve 
transits and e c lipses and the transit and eclipse from 
IKnutson et al.l (|2007() . we find that the scatter in the 
difference in centroids of the two stars (adding in quadra- 
ture the X and y components) ranges from 0.0035 to 
0.0042 pixels for the Gaussian centroid after binning by 
128 exposures (51 seconds); this is a factor of 3 — 7 times 
smaller than the flux-weighted centroid, and also no sys- 
tematic trend in x, and a weak systematic trend in y. 
Since the Gaussian centroids of the two stars track one 
another well and the difference in their positions is nearly 
uncorrelated with their position on the detector, we are 
confident that the Gaussian centroid is giving the cor- 
rect absolute position of these stars. When the data 
are fit with a 3.5 pixel radius aperture (as described 
in more detail below), the 2-D Gaussian centroid yields 
a x^ which is smaller for 9 of 14 transits/eclipses than 
the flux-weighted centroid, while the total x^ is smaller 
by 130 (after discarding the first 55 minutes of data for 
each eclipse/transit which has the steepest portion of the 
ramp). 

In conclusion, we recommend the 2-D Gaussian for 
Spitzer IRAC Channel 4 sub-array centroiding of bright 
targets as it appears to behave in a near-optimal manner. 
As we will show, this results in very small scatter in the 
resulting photometry. 

4. RAW PHOTOMETRY 

We carried out photometry on our data using aperture 
photometry with a circular aperture. The contribution of 
the pixels on the edges of the circle are calculated by mul- 
tiplying the total fiux in the pixel by the geometric frac- 
tion of the the pixel that is covered by the circular aper- 
ture. This is done using the GSFC Astronomy Library 
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Fig. 1. — Comparison of the centroiding algorithms for HD 
189733 versus M dwarf companion (top panels) and HD 80606 ver- 
sus HD 80607 (bottom panels) in the z-direction (left panels) and 
y direction (right panels). The black dots are for the flux- weighted 
centroid, the red dots for the 2-D Gaussian centroid, and the green 
dots for the parabolic centroid. The HD 189733 data have been 
binned by 512 exposures (205 seconds), while the HD 80606 data 
have been binned by 4 exposures (56 seconds). 



IDL routine pixwt.pro. Figure [2] shows a logarithmically- 
scaled median image from one of our sets of observations. 
As can be seen in the image, the sub-array is 32 pixels 
square, and a companion M-dwarf lies 9 pixels from the 
target star. We tried a range of apertures, discussed be- 
low in ij4.31 but our final analysis uses a 4.5 pixel radius 
aperture which is shown as a red circle centered on the 
target star. Note that this aperture size contains the 
bulk of the target flux, and is near the minimum in flux 
just inside the first Airy ring; this makes our photome- 
try less sensitive to variations in position. We have fit 
both stars with the measured point response function for 
IRAC Channel 4, and we find that the contribution of the 
M dwarf within this aperture is less than 0.06% of the 
target star flux for all of our observations. The resulting 
light curves for our twelve observatio ns plus the transit 
and eclipse from lKnutson et aD ()2007f ) are shown in Fig- 
ure [3] for an aperture of 4.5 pixels radius. Note that for 
each transit/eclipse pair, the flux at eclipse is higher than 
the flux at transit; this indicates the planet is brighter 
on the day side than the night side. 

4.1. Background subtraction 

To subtract the backgr ound, we used a similar proce- 
dure as that described in lKnutson et al.l (2007) : namely 
we fit a Gaussian to a histogram of the counts from a 
subset of pixels located in the corners of the image (ex- 
cluding the M dwarf companion, and excluding the top 
row). This background contributes about 1.9-2.6% of the 
total flux in our 4.5 pixel radius aperture, which we sub- 
tr act from the time s eries frame by frame. As discussed 
by [Harrington et al.l ()2007[ ) and in the IRAC Instrument 
Handbook, the flux and the background of the lst-5th 
and 58th frame of every set of 64 exposures is systemat- 
ically lower than the other exposures in a set. However, 
after we carried out background subtraction frame by 
frame, the offset in these exposures does not appear in 
our total time series. Consequently, we believe it is due 
to a bias offset that affects the entire frame uniformly, 
and thus is easily removed. 
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Fig. 2. — Logarithmic scaling of a medianed image from one of 
our observations; horizontal and vertical axes are pixels. White 
represents about 13,700 counts per pixel, while black is about 6 
counts per pixel. The red circle is a 4.5 pixel radius aperture. 



4.2. Detector ramp 

Spitzer was not envisioned as an instrument for carry- 
ing out high-precision photometry on bright targets; con- 
sequently it was designed without sub-mmag exoplanet 
photometry in mind. An instrumental artifact that ap- 
pears at the '^mmag level in photometry, but can be up 
to 10% for low-illuminatio n pixels over 3 3 hours, is the 
so-called "detector ramp" ([Demind 120091 ): a pixel which 
is illuminated uniformly in time shows a gradual increase 
in the detected flux (see Figure [3]). This is an important 
effect to correct for in fltting photometric time-series; 
unfortunately there has not been a full understanding of 
this effect for the Spitzer IRAC detector. Here we derive 
a toy-model which qualitatively matches the behavior of 
the ramp ( §4.2.ip . We show that prior functions used for 
ramp-corrections in other analyses of IRAC data (e.g. 
iDeming et alll2006t iDesert et al1l2009D do not have the 
correct functional form to describe the observed ramp 

(mm- 

Understanding how to remov e the detector ramp has 
evolved with time. Early work (jDeming et al.l[2C)06( ) ra- 
tioed the transit star to other sources, and modeled the 
baseline in the ratio as linear or quadratic in time. Fit- 
ting functions to t he ramp directly have u sed either ex- 
ponentials in time (iHarrington et al]|2007D or polynomi- 
nals in the log of time (jKnutson et al.ll2009l : IDesert et al.l 
|2009| ). These approaches have been adequate at lower 
signal-to-noise ratios, but for the present high S/N data 
we are motivated to flnd an improved functional form. 
We propose a new functional form, motivated by the toy 
model, which has the correct behavior and matches the 
observed ramps better (i i4.2.3|) . We apply a range of tests 
to this ramp model, and show that it performs better 
than other ramp models f §4.2.4|) . 
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6.50 b 



2000 



4000 6000 

Doto point (binned by 64) 



8000 



10000 



Fig. 3. — Atlas of transits and secondary eclipses obtained at 8 microns with Spitzer. The photon counts per exposure averaged over 64 
exposures is plotted versus the sequence of each set of 64 exposures. The numbers above or below each transit/eclipse indicate the orbital 
phase. The solid red curves show the best-fit double-exponential ramp model (which includes the phase-variation of the planet for the 
secondary eclipses). 



4.2.1. Toy model for the detector ramp 

The ramp effect is hypothesized to be due to trap- 
ping of electrons in detector defects ( "charge-trapping" ) . 
When a pixel is first illuminated, the charge traps are ef- 
fectively empty, and some fraction of the electrons gener- 
ated by the incident flux are retained by the traps instead 
of being read out by the array. As these charge traps fill, 
the effective gain of the detector goes up, until even- 
tually the effect disappears. Thus bright, non-variable 
sources should have a detected flux that asymptotes to 
a constant value. When a pixel is not illuminated (or 
illuminated at very low intensity), the trapped charge 
gradually releases with time, causing the charge traps 
to become empty; this leads to ghost images after ex- 
posure of bright sources. A consequence of this model 
is that higher illumination pixels fill their charge traps 
more quickly, thus showing a much shorter detector ramp 
timcscale. Although it is not clear that this model is 
correct, its predicted behavior agrees qualitatively with 
the observed IRAC photometric properties: for a bright 
source, the central pixels have a short ramp which satu- 
rates quickly, while the pixels with lower illumination in 
the wings of the PSF show a more gradual ramp. This 
behavior, though, is difficult to model quantitatively as 
the pixel illumination varies with time due to Spitzer 



pointing variations (see section [ 

A simple toy model can be developed for charge- 
trapping as follows. Let 7m be the fraction of volume 
of a detector pixel filled with charge-traps, j{t) be the 
fraction of volume of a pixel with empty charge-traps 
at time t, and /3 be the total well-depth (electrons). A 
pixel is illuminated below saturation with an intensity 
causing I{t) electrons to be released per second, while 
the measured intensity is I'{t) (e~ sec""'^). As the pixel 
is illuminated, the charge traps fill up at a rate propor- 
tional to the intensity times the fraction of volume of 
empty charge traps; however, there is also a timescale r 
at which electrons in charge traps are released, causing 
7 to increase. This gives 



dj _ _ljfl _ 7m - 

dt ~ (5 ^ T 

The measured intensity is then 



/'(t) = (l-7)/(t) + /3- 



7(0 



(1) 



(2) 



These equations have no closed-form solution for an ar- 
bitrary I(t); however, we can solve for their behavior in 
certain limits. For instance, if the illuminating intensity 
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is constant, I{t) = /q, for times t > to, then 

7(i) = 7m(l-/r//3)-i + (7(to) - 7m(l - It//3)) e(--^-V«(*- 

(3) 

As we are observing a bright star, we can simplify this 

equation by assuming that t~^ <^ -fo//3, but we are stiU 

below saturation and in the linear regime (/o^oxp ^ 0.9/3, 

where iexp is the exposure time). In this limit 

I'it) « /o(l - 7) - /o (l - 7(io)e"*'*~*°^) ■ (4) 

This gives the expected ramp behavior: more strongly 
illuminated pixels have an apparent intensity, /', that 
asymptotes to a constant more quickly on a timescale 
P/Iq. At modest illumination, this becomes 

I'{t) = /o(l - 70 + lohr'it - to)); (5) 

a gradual linear ramp. 
In the limit of zero illumination, 

7(i)=7„-(7„,-7o)e"(*"*°)/^ (6) 

Thus the apparent intensity is 

/'W-^(7m-7o)e-^*-*"^/^ (7) 

r 

This leads to persistent or ghost images that decay ex- 
ponentially with time after observation of a bright target 
when the illumination is so strong that 70 ^ 7m- 

4.2.2. Prior models for the detector ramp 

The correction for the detector ramp is typically ap- 
plied after performing photometry on the target star 
rather than at the pixel level, with some ex ceptions (e.g. 
iKnutson et all l2007t iLaughlin et"alll2009[ ). There is a 
simple reason for this: at the pixel level it is very difficult 
to disentangle the ramp from pointing variations, while 
aperture photometry with a sufficiently wide aperture 
gets rid of most of the pointing variations and isolates the 
ramp behavior. Most ramp corrections have simply been 
functions that match the behavior of the ramp; two popu- 
lar functions are ao-l-ai(t — io)+a2 log {t — to) (log-linear) 
and flo + fli log (i — io) + ^2 [log {t — to)] (quadratic log) 
whic h both seem to work well for IRAC Channel 4 data 
(e.g. iDeming et al.l 2006; Doming 2009). In particular, 
the logarithmic term matches the shape of the ramp well, 
which is steeper towards the beginning and shallower to- 
wards the end. 

Given the toy model described in the prior section, 
this logarithmic behavior would appear to be largely co- 
incidental. Aperture photometry combines pixels with 
a wide range of illuminations; those with high illumina- 
tion, which is most of the flux, have a ramp that is steep 
and saturates quickly, while those with low illuminations 
in the wings of the PSF have a longer timescale ramp 
that saturates more slowly. Summing up these short and 
long timescale ramps gives a shape which is steeper in 
the beginning and more gradual as time passes, which is 
well modeled by a logarithmic function. The linear term 
or a squared logarithmic term gives enough extra degrees 
of freedom to the model to adjust the slope of the curve 
and give a good fit to most ramp data. This model has 



the additional advantage that it is linear (except the ini- 
tial time, to, in the logarithmic term), and thus is quick 
*'^nd easy to fit to the observed ramp. 

However, the log-linear and quadratic log ramp models 
have a serious drawback: they do not have the correct 
behavior on long timescales. Both the log function and 
linear function increase without bound, while the detec- 
tor ramp does appear to saturate at a constant value for 
the brightest pixels. Thus, with a dataset with long du- 
ration, the log plus linear model or quadratic log ramp 
models should do a poor job in fitting the ramp shape. 
In addition, the log plus linear and quadratic log mod- 
els do not describe what the final asymptotic flux value 
will be, and thus does not give a ramp correction, but 
only gives an empirical flt to how the flux is varying 
over the timescale of a given observation. These points 
are particularly important for small aperture photome- 
try where most pixels have high illumination and thus 
saturate quickly. Consequently, we advocate not using 
ramps that are polynomials in time and/or log time. 

4.2.3. New model for the detector ramp 

Motivated by the toy model in section I4.2.1[ we de- 
cided to try an exponential ramp function. As this model 
predicts, the time constant is a function of pixel illumi- 
nation. However, due to the pointing variations, wc were 
not successful in correcting for the ramp on the pixel 
level. Instead we tried a ramp correction function that 
is simply the sum of two exponential terms: 



F'/F = ao- aie 



-t/n 



026 



-t/T2 



(8) 



where F' is the flux affected by the ramp, and F is the 
flux corrected for the ramp. Although this does not have 
exactly the correct behavior for the sum of pixels with 
different illuminations (assuming the toy model is cor- 
rect), it does have the correct asymptotic behavior, and 
qualitatively represents the correction from higher and 
lower illumination pixels. 

Figure [3] shows the ramp function overplotted with our 
data for the fourteen transits and eclipses. 

4.2.4. Performance of double- exponential ramp 

In addition to the qualitatively correct behavior of the 
double-exponential ramp, we flnd that this ramp function 
leads to a smaller scatter in our derived eclipse depths 
for the seven eclipses, as well as less sensitivity to the 
various choices we make in our analysis. We held the 
planet-star radius ratio and impact parameter flxed at 
the transit values when analyzing the secondary eclipses. 
We ran initial fits for each ramp on photometry com- 
puted for a 3.5 pixel radius aperture with the first 55 
minutes for each transit/eclipse discarded, and then de- 
termined how the eclipse depth changed as we varied 
individual analysis parameters. The scatter in the seven 
secondary eclipse depths for the double-exponential ramp 
model is smaller by 30% than for the log-linear ramp 
(3.05% versus 3.94%), and smaller by 20% compared to 
the quadratic log ramp (3.05% versus 3.68%). An addi- 
tional indication of the more robust behavior of this ramp 
function is that the mean depth only changes by 0.2% if 
we first fit the ramp to the out-of-transit/eclipse data, 
and then fit the transit/eclipse to the ramp-corrected 
data versus a simultaneous fit to the transit /eclipse and 
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ramp. The log-linear and quadratic log ramps have a 
mean eclipse depth that changes by 0.5% and 1%, respec- 
tively, between these two reduction techniques. Likewise, 
the double-exponential ramp changes in eclipse depth by 
only -1.1% if the first 55 minutes are discarded, while the 
log-linear and quadratic log change by -1.2% and 3.7%, 
respectively. The double-exponential ramp is also less 
sensitive to aperture size. For apertures between 3.5 to 
5.0 pixels in radius, the individual eclipse depths vary by 
1%, while for the log-linear and quadratic log ramps, the 
variation is 1.5% and 2.4%, respectively. Finally, the to- 
tal x^ for the double-exponential ramp model is slightly 
smaller by 21 than the log-linear ramp, and by 27 than 
the quadratic log ramp, which by the F-test for the ad- 
ditional 13 free parameters (for seven transits and six 
eclipses; the phase -function eclipse has no ramp) favors 
the double-exponential ramp at > 99.999% confidence. 

There are two drawbacks of the double-exponential 
ramp function: (1) it involves two non-linear fit param- 
eters, Ti and T2, which need to be optimized with a non- 
linear minimizer; (2) in some cases when there is very lit- 
tle ramp (possibly due to high illumination prior to our 
observations), one or both of the r values can diverge, 
or in some cases they can become degenerate. However, 
these drawbacks are straightforward to deal with by set- 
ting bounds in a non-linear solver, and are outweighed 
by the improved fit to the observed ramp, the correct 
asymptotic behavior, the smaller scatter in our results, 
weaker dependence on aperture size, the weaker depen- 
dence on whether the ramp is first corrected or fit si- 
multaneously, and less sensitivity to whether the steep 
portion of the ramp is discarded. Thus, we advocate 
using this ramp function for IRAC Channel 4 data. 

4.3. Aperture size 

We carried out photometry with apertures ranging in 
radius from 1 pixel to 7 pixels. We fit each transit and 
eclipse separately, and then computed the standard de- 
viation of the data divided by the best-fit model (this 
is essentially the residuals in magnitudes). We find that 
the residuals are minimized at 4.39 mmag for an aper- 
tu re radius of 3.5 pixels; this is the same aperture chosen 
in lKnutson et al.l (|2007[ ). For aperture radii of 3.0 pixels, 
the scatter is 4.5 mmag, while for 4.0 pixels it is 4.40 
mmag, and for 4.5 pixels is 4.47 mmag, indicating that 
there is a shallow dependence on scatter with aperture 
size. 

More importantly, we wish to minimize the presence of 
red noise in the residuals of the data. Consequently, we 
looked at the scatter in the residuals binned over a range 
of bin sizes from one exposure to 1920 exposures as a 
function of aperture size. We then took the mean of the 
scatter of the binned data over all fourteen transits and 
eclipses, and computed the product of this mean scatter 
divided by the unbinned scatter over all bin sizes. The 
minimum occurs for an aperture of 4.5 pixels; although 
this has a slightly larger residual scatter without binning, 
the binned residuals are smaller relative to the unbinned 
residuals than for the 3.5 pixel radius aperture case. We 
have also measured the power spectrum of the residuals, 
and we find that the 4.5 pixel radius aperture minimizes 
the long period power. Thus, we feel that this aperture 
size represents an appropriate compromise between small 
scatter in the unbinned data (which varies weakly with 



aperture size) versus minimization of the red noise com- 
ponent. 

Figure |4] shows the scatter in the binned residuals, av- 
eraged over the 14 transits, as a function of bin size. Even 
up to bin sizes of 8832 exposures (1 hour bin), the scatter 
in the data does not deviate significantly from the inverse 
square root of the bin-size; this indicates that the resid- 
uals are uncorrected, and thus there is little (if any) red 
noise present. Remarkably the scatter in the one hour 
bins reaches 30 /xmag; however this is after subtracting 
off the double-exponential ramp model. 




Fig. 4. — Scatter in the data (vertical axis) divided by the model 
binned by a number of bins (horizontal axis), averaged over the 14 
transits and eclipses presented here. Red line is the extrapolation 
from the unbinned data by the inverse square root of the bin size. 



For an aperture radius of 4.5 pixels, the median counts 
per exposure is 66,792. If photon counting errors domi- 
nate, then the expected noise level is 3.87 mmag per ex- 
posure. Including read noise (4.5 e~ per pixel) and sky 
noise (^30 counts per pixel), the expected uncertainty is 
3.99 mmag (we did not use the BCD uncertainties since 
these overpredict the noise properties according to the 
Spitzer Observer's Manual). The standard deviation in 
the residuals is 4.47 mmag, which is only 15.5% greater 
than the photon counting error and 12.0% greater than 
the expected errors including read noise and sky noise. 
Thus, the noise properties after correcting for the de- 
tector ramp are very close to the expected photon noise. 
The 3.5 pixel radius aperture has a residual scatter which 
is only 11% above the photon noise; however this aper- 
ture size appeared to have more significant red noise, so 
we opted for the larger 4.5 pixel radius aperture. 

4.4. Conversion to barycentric Julian date 

We use the JPL Horizons ephemeris for the Spitzer or- 
bit to convert the satellite time (keyword DATE_OBS) 
to Barycentric Julian Date (BJD) in Coordinated Uni- 
versal Time (UTC).^ This correction is important since 
heliocentric and barycentric Julian Date can differ by as 
much as a few seconds, and different time systems can 
vary by seconds depending on the number of leap sec- 
onds included, which is close to the timing precis ion we 
can achieve with these data (jEastman et al.ll2010[ ). 

* There is an additional -1-65.184 sec offset to convert to Barycen- 
tric Dynamical Time for our data. 
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4.5. Error analysis 

We compute the errors on model parameters by calcu- 
lating the residuals from each fit, shifting these by a ran- 
dom number of observations, and then adding the shifted 
residuals back to the best-fit model, t hen re-fitting; the 
so ca lled "prayer-bead" analysis (e.g. lAgol and SteffenI 
l2007t ). For each transit and eclipse we carried out 2000 
iterations of the prayer-bead shifts. This has the advan- 
tage of preserving correlations in the noise of the data 
that might still be present. For instance, if the ramp- 
model is incorrect, there may be systematic deviation due 
to using the wrong ramp model, and these deviations are 
preserved within the residuals. This approach has some 
disadvantages; for instance, if the noise behaves differ- 
ently within eclipse than outside eclipse, this might ex- 
aggerate the noise outside eclipse. Another disadvantage 
of this technique is that the number of independent trials 
is limited by the size of the data set since point-order has 
to be preserved; consequently we also randomly chose to 
reverse the residuals or change their sign to give more 
independent noise realizations. There is also the possi- 
bility that the effects of correlated noise may be removed 
in the fit. Even with these disadvantages we expect that 
this technique gives a fairly conservative estimate for the 
uncertainties on model parameters. 

5. FIT FOR STELLAR AND PLANET VARIABILITY 

Stellar variability can affect our fits to the transits and 
eclipses, as well as our estimate of the planet y a riabil - 
ity. We follow t he ap proaches in iKnutson et al.l ()2009f ) 
and ISing et all ()2009D to derive a new estimate of the 
relation between the optical and 8 micron flux variabil- 
ity of_the_stai__b^ comparing our data set with that 
of iHenrv and Wimi (2008), plus additional unpublished 
data. The contemporaneous optical monitoring data 
were taken with the TIO 0.8 m automated photometric 
telescope at Fairborn Observatory, which has a median 
time sampling of 1 day, but gaps of up to 2 weeks. The 
optical time series consists of the mean of Stromgren b 
and y magnitudes, subtracted from a nearby comparison 
star, HD 189410, giving /opt = A[(6 + y)/2 ] as a func- 
tion of time, as described in more detail in (jWinn et al.l 
l2007t ). Data outliers are removed (usually taken in poor 
conditions) , resulting in a total of 700 observations over 
5 seasons. 

Using the measured s tellar rotation pe r iod o f P* = 
11.953 ±0.009 days from iHenrv and Wind (pOOl . we fit 
a sinusoidal function to data within each season to in- 
terpolate the measured optical fluxes to the times of our 
Spitzer observations. For all but one of our IRAC obser- 
vations there were at least one optical observation taken 
within 1 day, and all within 2 days. We then computed 
the total unocculted 8 micron flux, /ir, at the mid-transit 
and mid-eclipse times of our 14 observations, after cor- 
recting for the best-fit double-exponential ramp, to look 
for a correlation between the 8 micron and optical fluxes. 
The initial data seemed to show little correlation between 
the infrared and optical fluxes, so we carried out a regres- 
sion of the infrared fluxes against flve variables: (1) the 
optical flux, /opt, which is entirely due to the star; (2) 
the phase $ which determines whether the source is tran- 
siting or eclipsing (i.e. whether we are seeing the day or 
night side of the planet), $ = during transit and $ = 1 



during secondary eclipse; (3) the average centroid posi- 
tion, Xc, on the detector for each of our observations (the 
y position varied little between observations); (4) the av- 
erage infrared background flux scaled to our aperture, 
/bkd (this was already subtracted in earlier analysis of 
the data, but we nevertheless include it in the regression); 
and (5) the amplitude of the first exponential ramp, oi. 
The last three terms are included to take into account 
the possibility of flat-field errors, imperfect background 
subtraction, and the imperfect performance of our ramp 
function. 
We find the best-fit relation 



JlL. 

(/ir) 



-1 = (0.197± 0.022) 



/opt \/opt/ 
(/opt) 



(1.044 ± 0.026)y^ - (6.59 ± 0.. 

(1.19 ±0.16) X 10-3* 
(1.55 ±0.48) X lO^'^ai 



X lO^^Xc 



(12.7 ±1.2) X 10" 



(9) 

where (/h) is the ramp-corrected flux averaged over all 14 
observations, while (/opt) is the average over the optical 
flux at the times of the 14 observations. The left hand 
side of this relation is plotted versus the right hand side 
in Figure[S] the scatter about this relation is 0.35 mmag. 
We have computed the uncertainties on the regression 
coefflcients by Monte Carlo simulation. 

The standard deviation of the residuals of our sinu- 
soidal fits to the optical data is 2.5 mmag (after exclu- 
sion of a few outliers), w hich is 1.8 times the opt ical flux 
uncertainty (1.4 mmag: iHenrv and Winnll2008[ ). Using 
the sinusoidal flts, the uncertainty on the optical flux we 
predict at the mid-points of our observed transits and 
eclipses ranges from 0.6-1.4 mmag after inflating the op- 
tical errors by a factor of 1.8. Since infrared stellar fluc- 
tuations are 20% of the optical, this predicts a scatter 
of 0.1-0.3 mmag in the infrared, which is consistent with 
the measured scatter. 

We have computed the expected spectral change at 8 
microns compared to (5 + y)/2 for a star spot model 
in which the star spots ar e modeled as Kurucz stellar 
atmospheres (|Kurucdll99l ) with 4000-4500 K (which is 
the estimated temp erature frorn occu lted star spots mea- 
sured with HST bv lPont et alll2008[ ). while the bulk of 
the st ar is 500 K, fo llowing the procedure described in 
Knuts on et al.l (|2009( ). We flnd the expected change at 8 
microns is 21-23% of the change in the optical, very close 
to our measured value, the flrst coefficient in equation [HI 

There are several important implications of this re- 
lation: (1) as the scatter in this relation is only 0.35 
mmag, this indicates that photometry with Spitzer is re- 
producible to 0.35 mmag over a 590 day period; (2) this 
scatter limits our uncertainty on the measurement of the 
night-side planet flux (during transit), so we can claim 
that the night side variability is less than 0.35 mmag, 
which is about 17% of the planet's night-side flux or 10% 
of the planet's day-side flux (which we flt for in section 
16. 7p : (3) the transits are 1.19±0.16 mmag fainter than the 
eclipses — after accounting for stellar variability — due 
to the cooler night side of the planet tha n the day side, 
confirming the phase variation detected in lKnutson et al.l 
(|2007f ): (4) the infrared flux variations are about 20% of 
the optical variations. 
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Fig. 5. — Left hand side of equation |9] versus the right hand side. 
Solid line is equality; the vertical scatter in the residuals of this 
relation is 0.35 mmag. 

The derived infrar ed/optical corre lation is nearly twice 
the value derived in lKnutson et alj (|2009), which used a 
smaller subset of data to carry out the correlation and 
thus was unable to regress against these other factors. 
Our estimate of the expected flux variations from Ku- 
rucz stellar atmospheres i ndicates that ou r deriv ed value 
is likely correct. However. iKnutson et al.l ()2009D derived 
a larger change in the stellar flux — by interpolating the 
observed y-band light curve — than we obtained by si- 
nusoidal fitting of the (6 + y)/2 light curve over the pe- 
riod of duration of the phase function observation, so the 
resulting estimates o f 8 micron stellar variation for the 
IKnutson et al.l ()2007() observation are nearly the same: 
a 0.6 mmag increase in stellar flux between transit and 
eclipse. 

6. ECLIPSE AND TRANSIT MODELS 

We fit a model of a straight-lined trajectory of the 
planet over the disk of the star. To compute the tran- 
sits and eclipses, we used the analytic formulae from 
iMandel and Agol ([2002), treating the planet as a uni- 
form disk (no limb-darkening), and the star as a disk 
with a linear limb-darkening law. 

For each transit the model has six physical parame- 
ters and four ramp parameters: the stellar flux F^:] the 
sky velocity v (units of stellar radius per day); the im- 
pact parameter b (units of stellar radius); the planet- 
star radius ratio p = Rp/R^ (dimensionless) ; the time 
of central transit t^ (in Barycentric Julian Day, BJD); 
the linear limb- darkening parameter of the star ui (di- 
mensionless); and the double-exponential ramp param- 
eters ai , Ti , 02 and T2 (equation [8]) . Note that we ne- 
glect the contribution of the planet's flux during transit; 
this is because wc find this is completely degenerate with 
the transit depth, and thu s leads to problems in fitting 
(jKipping and Tinettill2009() . We initially neglected phase 
variation of the planet and variation of the fiux of the star 
during transits as the ramp affects all of the transit data 
sets and thus a short timescale (5 hour) planet or stellar 
variation can not be disentangled from the ramp for a 
single observation. 

For the secondary eclipses, we assumed that the 
pl anet phase-function f ollowed the same shape as that 
of IKnutson et al.l (|2007D . which we held fixed in our fits 
to each secondary eclipse, but we allowed the total planet 
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Fig. 6. — (a) Average light curve for the seven transits and (b) 
seven eclipses with best- fit models (solid lines). The data have 
been binned by 8960 data points to a total of 69 data points for 
each. 

flux to vary for each eclipse. For each secondary eclipse 
we held fixed the planet/star radius ratio p, the im- 
pact parameter 5, and the velocity v, at the best-fit val- 
ues from the transit observations; these parameters are 
poorly constrained by the secondary eclipses, and hold- 
ing them fixed has no impact on the fitting. Thus, for 
each secondary eclipse we have three physical parameters 
that are varied: the stellar flux F^,; the planet fiux Fp-, 
the central time of eclipse tc] as well as the four ramp 
parameters. 

6.1. Results from fits to individual transits/eclipses 

We allowed the model parameters to vary indepen- 
dently for each transit/eclipse. These fits were necessary 
since a simultaneous fit to the entire data set is com- 
putationally intensive due to the large number of data 
points; we avoided pre-binning the data to preserve as 
much information as possible about the noise in the final 
results. Figure |6] shows the average of all seven transits 
and all seven eclipses, corrected for the detector ramp 
and folded to the same orbital phase. We have binned 
the data by 8960 exposures to 69 data points for clarity. 

6.2. Transit impact parameter and sky velocity 

For the transits, we found that the sky velocity, v, 
and impact parameter, 6, have no evidence for varia- 
tion. Figure [7] shows each of these parameters plot- 
ted versus the transit number. The sky velocity has 
an average fractional uncertainty of 0.62%; the scatter 
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Fig. 7. — (a) Impact parameter and (b) transit velocity, versus 
transit number with estimated error bars. The horizontal solid line 
is the average of the seven measurements, while the dotted lines 
are the uncertainties on the average values. The numbers above 
each data point show the number of periods before or after the 
zero-point of our measured ephemeris. 

in the measured values is 0.57%, and thus is consistent 
with being constant. Combining our data together, we 
find the average sky velocity is 25.125±0.064 i?* day~^, 
while a limit on the variation of the sky velocity is 
dv/dt = (-5.5 ± 6.6) x lO"'*^?* day~^. The impact pa- 
rameter has an average measured value of 0.6631±0.0023 
i?*, with an average fractional uncertainty for each ob- 
servation of 0.93% and a fractional scatter for the seven 
observations of 0.67%; also consistent with being con- 
stant. We constrain the change in impact parameter to 
be db/dt = (-0.02 ± 2.67) x lO'^R^ day-^. Thus our 
data indicate that the impact parameter and sky velocity 
of the transits remain constant to <1% over a duration 
of 590 days. 

6.3. Transit and eclipse times 

We measured the transit and eclipse times for the seven 
transits and eclipses, shown in Figure [8l as well as in 
Tables ll|2l The errors on the transit times range from 
2.4-3.6 seconds and are some of the most precise transit 
times ever measured; comparable to, or better than, th e 
three HST transit times reported in iPont et "all ()2007D . 
We fit separate ephemerides to the transits and eclipses; 
the results are shown in Table [3] If we instead fit the 
transit times with the quadratic function: tn — to + Pn + 
^PPn^, where tn is the time of the nth transit, and P is 
the change in period of the orbit, we find P = —0.06 ± 



0.02 sec yr~^. Since this is primarily due to the last data 
point, which may be an outlier, we do not view this as a 
significant detection. 

The uncertainty on the transit times and eclipse times, 
as well as the derived ephemerides, are inversely propor- 
tional to depth of the transits and eclipses. This is due to 
the fact that the timing precision is proportional to the 
inverse of the flux gradient with time during ingress and 
egress. The ingress and egress duration are the same 
for the transit and eclipse (assuming a circular orbit), 
so the ratio of the fiux gradient scales with the ratio of 
their depths. The ratio of the depths is proportional to 
the ratio of the surface brightness of the star to the sur- 
face brightness of the planet (limb-darkening is weak for 
this star at 8 microns), so the transit time precisions are 
smaller than the eclipse precisions by the ratio of the sur- 
face brightness of the planet to the star, which is about 
14.3%, or a factor of 7.0. 

Figure [5] shows the deviations from the ephemeris for 
both the transits and the eclipses. The transits have a 
scatter of 5.1 seconds, which is very close the the observa- 
tional errors; there is no evidence in our data for transit 
timing variations over a period of 590 days. The eclipses 
also appear to be precisely periodic - their scatter with 
respect to the best-fit ephemeris is 27 seconds, which is 
comparable to the errors on each data point. The period 
derived from the transits differs from that derived from 
the eclipses by only 0.1 seconds! 

The eclipses appear 69±11 seconds later than 1/2 of 
an orbital per i od aft er the transits. As discussed in 
iKnutson et al.l ()2007f ), this is partly due to the light- 
travel across the system, 30.8±0.6 seconds (this uncer- 
tainty is due to t he uncertainty in stellar mass, M* = 
0.806 ± 0.048Mr7.. [Torres et al.lf2008l ). while the remain- 
ing 38±11 seconds can be mostly accounted for by the 
hot spot on the planet causing an offset in the time of 
eclipse when the planet is modeled as a uniform disk, as 
shown in section [ 



6.4. Limits on the presence of companion planets from 
transit-timing 

These transit data show no significant timing varia- 
tions, but from these we can constrain the maximum 
mass allowed of additional planets in the system. Transit 
timings are a particularly sensitive probe for planets in or 
near mean-motion resonance (MMR) and previous stud- 
ies have ruled out Earth mass or super-Earth mass plan- 
ets in low-order MMR for several systems. Prior transit 
timi ng variation (TTV) analyses of the HD18973 3 sys- 
tem (jHrudkova et al.ll2010HMiller-Ricci et al.ll2008[ ) used 
data with timing precision of order 30 seconds and were 
sensitive to planetary masses of near (and below) one 
Earth mass in favorable MMRs. Our Spitzer observa- 
tions of HD189733 have nearly a factor of 10 better tim- 
ing precision and consequently have improved sensitivity 
to secondary planets by that same factor. Here we calcu- 
late the maximum mass that an additional planet could 
have based upon these transit data. To do so we note 
that the x^ P^r degree of freedom of the timing residuals 
is slightly more than unity. We therefore scale the timing 
uncertainty by a factor of 1.5 and then multiply by two 
to achieve our 2a (~ 95% confidence level) upper bound 
on the timing variations. 

Figure [9] shows 95% confidence-level constraints on sec- 
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TABLE 1 

Transit parameters 



Phase F* tc at^ Ate b v {RpjR,Y 
(counts) (BJD-2454000 days) (sec) (sec) (lU) (R* day~^) (%) 



"1 



— ai, —a2 
(counts) 



(10- 



,T2 
^days) 



-109.0 67222.6 37.611919±0.000034 3.0 -4.2 0.6694±0.0062 25.02±0.15 2.4088±0.0037 

1.0 66782.7 281. 655291±0. 000040 3.4 -0.0 0.6586±0.0066 25.28±0.16 2.4022±0.0047 

2.0 66722.4 283. 873884±0. 000042 3.6 1.4 0.6592±0.0066 25.24±0.18 2.4253±0.0063 

52.0 66758.8 394.802711±0. 000028 2.4 5.2 0.6636±0.0067 25.05±0.17 2.4333±0.0051 

63.0 66995.3 419. 207003±0. 000036 3.1 1.7 0.6594±0.0049 25.18±0.11 2.4225±0.0049 

158.0 67385.6 629. 971694±0. 000033 2.9 1.8 0.6641±0.0062 24.90±0.16 2.3984±0.0062 

159.0 67242.2 632. 190128±0. 000039 3.4 -10.4 0.6685±0.0060 24.99±0.16 2.3965±0.0074 



0.08± 0.03 
O.llit 0.03 
O.llit 0.03 
0.13± 0.03 
0.12± 0.02 
0.16± 0.02 
0.08± 0.03 



75, 550 
281, 107 
756, 752 
756, 712 
773, 722 
639, 570 
715, 669 



0.6010, 
0.4735, 
0.4929, 
0.4891, 
0.5322, 
0.4560, 
0.4890, 



6616 
5.4525 
5.5313 
5.6423 

2047 
5.5904 
5.8697 













TABLE 2 
















Eclipse parameters 




Phase 


F* 


tc 


Ttc 


At, 


Fp/F, 


— ai, — a2 


n,T2 




(counts) 


(BJD-2454000 days) 


(sec) 


(sec) 


(%) 


(counts) 


(10- 2 days) 


-108.5 


67466.9 


38.722278±0.000265 


22.9 


9.1 


0.3345±0.0057 


0, 


0.0000,0.0000 


0.5 


66870.3 


280. 546423±0. 000263 


22.7 


-32.5 


0.3469±0.0060 


609, 535 


0.3793,4.5632 


1.5 


66808.5 


282. 765713±0. 000350 


30.2 


29.3 


0.3420±0.0068 


255, 118 


0.1134,4.8741 


51.5 


66904.4 


393. 693798±0. 000414 


35.8 


-26.2 


0.3368±0.0042 


750, 727 


0.6350,6.1896 


63.5 


67097.2 


420. 317184±0. 000287 


24.8 


16.2 


0.3623±0.0064 


759, 647 


0.5089,6.0830 


157.5 


67582.4 


628. 862649±0. 000390 


33.7 


-30.8 


0.3378±0.0091 


330, 239 


0.0310,5.6520 


158.5 


67318.3 


631. 081875±0. 000313 


27.0 


25.5 


0.3477±0.0075 


737, 685 


0.5505,5.6994 



TABLE 3 
Transit and eclipse ephemerides 



To 
(BJDuTc) 



P 

(days) 



Transit: 
EcHpse: 



2454279.436714±0.000015 
2454279.437510±0.000125 



2. 21857567±0. 00000015 
2. 21857456±0. 00000131 



ondary planets with near circular orbits in this system 
based upon these data and the radial velocity (RV) mea- 
surements from Boisse 2009. T hese limits are de rived 
from the analytic formu la given in lAgol et al.l ()2005[ ) and 
iSteffen and Agoll (|2005r). We do not attempt an in-depth 
numerical analyis of these transit times here — the robust- 
ness of li mits derived from analyti c formulae was de mon- 
strated inlA gol and Stefrenl|2007D.| Nesvornvl (|2009f ) and 
iNesvornJa nd Bcaugc (201(|). These data exclude plan- 
ets above 2 Earth masses for any orbit that lies closer 
to the known planet than either the interior or exterior 
2:1 MMR. The transit-timing mass exclusion is supe- 
rior than the exclusion from radial- velocity data for pe- 
riods from 1 to 5 days, excluding all planets with masses 
greater than 3 Earth masses within this range. In ad- 
dition they exclude planets with masses well below the 
mass of Mars — approximately 0.2 Mars masses or 2 Moon 
masses at 95% percent confidence — in circularly orbiting 
2:1 or 3:2 MMRs (interior or exterior). For non-circular 
orbits the sensitivity generally increases. However, in 
low-order MMR the mass sensitivity can decline as much 
as a fa ctor of 10 for eccentric orbits — (see, for example, 
lAgol a nd Stcffcn 2007). 

6.5. Effect of hot spot on secondary eclipse time 

As discussed in iKnutson et al.l ()2007[ ). the 8 micron 
phase function indicates that the hottest point on the 
planet is offset from the sub-stel lar point. This was pre- 
dicted bv lCooper and ShowmanI (2005), and is attributed 
to the advection of energy by a super-rotating wind en- 
circling the equator of the planet. This offset hot spot 
means that the ingress and egress of the secondary eclipse 
will have a shape that differs from our model, which 



utilizes a uniform disk. In particular, this means that 
the steepest portion of ingress and egress will be offset 
from the uniform disk case; since the hotspot is on the 
trailing side of the planet with respect to the direction 
of motion, this causes a delay in the eclipse time when 
fit with a uniform di sk m odel (iCharbonneau et aLll2005l : 
I Williams etHI 120061 ). In IKnutson et all (|2007D we esti- 
mated that the hot spot would cause a delay of at most 
20 seconds; however, the fit to the phase function in that 
paper did not correct for stellar variation which caused 
the location of the hot spot to be underestimated, leading 
to an underestimated uniform time offset. 

To estimate the magnitude of this effect, we used a sim- 
plified model of the longitudinal p lanet b rightness which 
is discussed in iCowan and Agoll (|2010d ). Briefly, each 
position on the planet is treated as a parcel of gas which 
moves eastward, absorbing star light as it passes across 
the day side, all the while radiating with a time constant 
Trad- This model can be parameterized by a single pa- 
rameter, e = Trad/Tadv, where Tadv IS the time it takes 
a parcel of gas to circle the planet. Small values of e 
("instant" re-radiation) lead to darker night sides and 
day side temperatures which are in equilibrium with the 
incident stellar flux. Large values of e lead to nearly uni- 
form temperatures at each latitude. Thus in the small 
or large e limits we expect no timing offset since the day 
sides are symmetric. Only with e ~ 1 is there an off- 
set hot spot causing a phase function which peaks before 
secondary eclipse, as well as a slight offset in the times 
of eclipse ingress and egress if fit with a uniform planet. 

We computed the effect of e on the time of secondary 
eclipse by solving for the planet day-side longitudinal 
surface brightness at the equator in the Rayleigh- Jeans 
limit and assuming a constant temperature with latitude. 
We computed the eclipse ingress and egress from this 
model for the planet surface brightness, we fit this sim- 
ulated eclipse light curve with a model for the eclipse of 
a uniform planet, and from this best fit we determined 
the offset in the time of eclipse, the so-called "uniform 



12 



Agol et al. 




0.001 F 



Transit number 




Transit/Eclipse 



Fig. 8. — (a) Transit timing variations, observed minus cal- 
culated for a constant ephemeris, O — C, and (b) both transit 
and eclipse timing variations (ETV), O — C, versus transit/eclipse 
chronological order with estimated error bars. Note that panel (b) 
has a vertical scale that is 12 times larger than panel (a). The 
horizontal dotted lines are the average of the seven transits and 
seven eclipses; this is zero for the transits as we have subtracted 
off the best-fit transit ephemeris. The numbers above each data 
point show the number of periods before or after the zero-point of 
our measured ephemeris; the points are not plotted as they occur 
in time, but are simply evenly spaced. 



time offset" defined by IWilliams et alj (|2006[) . Figure 
[TUl shows this time offset as a function e; a positive off- 
set means that the secondary echpse occurs later than 
expected for a uniform planet. The maximum offset pre- 
dicted by this model is 43 seconds, which agrees with 
the observed eclipse time offset. Our measured eclipse 
time is plotted as a dashed line in this figure, with the 
uncertainty indicated by the horizontal shaded rectangle. 
Th e location of th e peak in the planet phase function 
from ' Knutson et al.l ()2007[ ) provides another constraint 
on the location of this hot spot, or equivalently on the 
value of e (see Figure [TT]). We used the relation between 
the infrared and optical stellar variability derived in sec- 
tion [5] to derive the change in stellar flux at 8 micron 
during the phase function measurement, about 0.6 mmag 
over 26.6 hours. We then fit the last 2/3 of the measured 
8 micron phase function to estimate e (Figure fTTI) : we dis- 
carded the first 1/3 of the phase function data since it is 
strongly affected by the ramp correction. We find a best- 
fit value of e = 0.74 ± 0.07, which we have also plotted 
as a vertical shaded region in Figure [TUl This value of e 
predicts a timing offset of 33 seconds, which is consistent 
with the measured offset of 38 ± 11 seconds. Figure [12] 
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Fig. 9. — Constraints (95% confidence level) on initially circular 
orbiting secondary planets in HD189733 as a function of the pe- 
riod ratio of the known planet based upon these transit data and 
the radial velocity measurements presented in Boisse 2009. The 
dotted curve are the limi ts from a TTV a nalysis alone from equa- 
tions (A7) and (AS) in .Agol et al.1 II2005I ). The dashed line is the 
expected sensitivity from 33 RV measurements with 3.5 m/s preci- 
sion calculated using equation (2) from Steffen & Agol (2005). The 
solid curve is the overall sensitivity from both RV and TTV mea- 
surements (summed in quadrature); the region above this curve is 
excluded. The solid dots are the variation in mean-motion reso- 
nance, ^ PtTansmp/mt, where mt.p are th e masses of the transiting 
and perturbing planets IIAgol et al.ll2005l ). Finally, the horizontal 
dot-dashed and triple-dot-dashed lines correspond to the mass of 
the Earth and the mass of Mars, respectively. 

shows a direct comparison of the secondary eclipse to the 
average of our seven secondary eclipses. The top panel 
shows the binned data as well as the best-fit secondary 
eclipse at 1/2 orbital period after transit plus the 30.8 
second light travel time delay (solid line), as well as the 
e = 0.74 model with light-travel time delay (dashed line). 
The bottom panel shows the residuals binned into eleven 
bins: pre-ingress, post-egress, eclipse, and four bins each 
in ingress and egress; the residuals are plotted for the 
uniform planet model (diamonds) and e — 0.74 model 
(filled circles with error bars). The uniform planet model 
shows points which arc on average higher in ingress and 
lower in egress, which is a sign of the shifted hot spot. 
The hot spot model provides a better fit to the data, al- 
though there is still scatter in the residuals which just 
reflects the low significance of the eclipse hot spot detec- 
tion (the uniform time offset is only 3.5(7: 38± 11 sec). 
Note that we have not optimized the hot spot model, 
but only computed the light curve from the best fit to 
the phase function (which gives e = 0.74). 

Consequently, there is no evidence for non-zero e cos w 
in this system. For the estimated value of e, the remain- 
ing time offset is 6 ± 11 seconds, which yields ecosw = 
0.00005 ± 0.00009. 

Another prediction of this model is the day-night 
brightness difference. For e = 0.74, we predict a night- 
side brightness which is 57% of the day-side brightness 
(day defined at mid-eclipse; night at mid-transit). This 
corresponds to a decrease in brightness from the day to 
night side which is 0.15% of the stellar fiux, which is 
very close to the value of 1.2±0.2 mmag derived in sec- 
tion [5] The minimum planet brightness (for the visible 
hemisphere) divided by the maximum planet brightness 
for this best- fit model is 50%, while the peak in observed 
planet brightness is 23 degrees, 0.065 orbital periods, or 
3.5 hours before the secondary eclipse. On the planet. 
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the hottest longitude is 13 degrees from the sub-stellar 
point; note that this differs from the hemispherically in- 
tegrated peak due to asymmetry in the longitudinal day- 
side intensity. Figure [TT] shows the phase function data 
after correction for stellar variation and binned by 1000 
data points, versus planet orbital phase with the best-fit 
model for the planet variability, e = 0.74, overplotted. 
Also plotted is our estimate of the night side flux based 
on equation [9] 




10.0 



Fig. 10. — Timing offset for a hot spot model as a function of the 
ratio of the radiative to advective time scales. The dashed line is 
the best-fit eclipse time offset after correction for light travel time, 
and horizontal rectangular shaded region is the l-cr confidence limit 
on this time. The vertical rectangular shaded region is the best-fit 
value of e to the 8 micron phase function, after correcting for stellar 
variability. 
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Fig. 11. — Planet phase function after correction for stellar vari- 
ability versus planet orbital phase. We only use the last ~ 2/3 
of the phase function to avoid problems with the ramp correction, 
and we masked the secondary eclipse. The thick solid line is the 
best-fit model for planet variability with e = 0.74. The dot with 
error bar on the left is our estimate of the night-side brightness 
from equation |9] The dotted line shows our correction for stellar 
variability during the phase function. 



6.6. Transit depth variation 
In the fits to the transits we allowed the depth of each 




Fig. 12. — Plot of average of 7 eclipses (top panel) with best-fit 
uniform planet model, offset by 30.8 seconds after 1/2 orbital pe- 
riod after the transit ephemeris (solid line), as well as the e = 0.74 
model with light-travel time delay (dashed line). The bottom 
panel shows the residuals binned into eleven bins: pre-ingress, 
post-egress, eclipse, and four bins each in ingress and egress; the 
residuals are plotted for the uniform planet model (diamonds) and 
e = 0.74 model (filled circles with error bars). 



planet to stellar radius, p = Rp/R*. This ratio also af- 
fects the duration of ingress and egress, but the amount 
of time spent in ingress or egress is small, so the primary 
effect is on the depth of transit. Figure [T51 shows the de- 
rived transit depths, p^, for the seven observed transits. 
There is evidence for variability in the transit depth - the 
uncertainty on the individual transit depths ranges from 
37-74 /xmag, while the scatter is 145 /zmag. This corre- 
sponds to a scatter in the fractional variation of transit 
depth of 0.6%, while the ratio of maximum to minimum 
is 1.5%. Fitting the transit depths with a single value 
of p gives a Ax^ of 40.8 for 6 degrees of freedom. Thus, 
transit depth variation is detected with high significance. 

We allowed the limb- darkening parameter to vary for 
each transit, which ranged from 0.08-0.13 for the seven 
transits, with an overall mean of ui = 0.12 ± 0.01. 
Even though limb-darkening is wea ker in the inf rared, 
an LTE Kurucz model atmosphere (jKurucd 11992') with 
parameters close to the values inferred for HD 189733, 
T^ff = 5000 K, [Fe/H] = 0.0, and I o g[g {cm /s^)] = 4.5, 
predicts a linear limb-darkening coefficient of 0.136, close 
to what we measure. We checked that the variation in 
transit depth is not due to variations in the best-fit limb- 
darkening parameter by holding the limb-darkening co- 
efficient fixed at the mean value; this did not affect our 
measured values of transit depth. 

In fact, the variation in transit depth is not necessarily 
due to a change in p. The average depth of each transit 
is giv en by S = (/path)7ri?p/(F* + Fp) (jMand el and Agoll 
120021 ) ■ where (/path) is the average surface brightness 
within the path of the planet across the star and Fp, F^ 
are the planet and stellar flux during transit. Although 
our model assumes a linear limb-darkening law, if the 
path of the planet passes over a region of the star with 
brighter than average surface brightness, then a larger 
depth will be inferred. 

So, it is possible that the change in transit depth is due 



to changes in (/path)//"*, Rp, R*, or Fp/F^. Variations in 
transit to vary independently through the ratio of the Rp or R^ seem unlikely to be responsible for the transit 
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depth variation as this would require changes in radius 
of 0.3% for either body: either 220 km for the planet {^ 
1/3 of a thermal scale height) or 1600 km for the star. 
Both of these variations seem too large to be physically 
plausible. Variations in (/path)/^* due to fluctuations in 
F^ are ruled out as the transit depth variations are un- 
correlated with the optical stellar flux variations. Vari- 
ations in Fp (the night-side planet flux) are less than 
0.35 mmag relative to i*"*, which is too small by a factor 
of 17 to be responsible for the transit depth variations. 
Thus, the most likely possibility is that the transit depth 
variations are due to a variation in the occulted stellar 
intensity, (/path)- This requires only a variation of 0.6% 
in the surface brightness of the path of the planet relative 
to the average stellar surface brightness, which is much 
smaller than the 12% change in surface brightness from 
center to limb inferred for the best-fit limb-darkening. 
We have checked that individual star spots are not re- 
sponsible for this variation by computing the standard 
deviation of the residuals in transit divided by the square 
root of the counts, which is 1.153, versus the same quan- 
tity computed out of transit, which is 1.154, so there is 
no evidence for individual star spots causing this differ- 
ence. Si milar inferences have b een drawn for Channels 1 
and 3 bv lBeaulieu et al.l (|2008[ ). while star spots are eas- 
ier to detect in the optical where the contrast between 
star spots and the stellar dis k is larger, and ha ve already 
been detected for this star bv lPont et al.l (|2007f ). who also 
resolved the shape and color-dependence of the spots. 




Fig. 13. — Transit depths measured for seven transits. Horizontal 
solid line measures the weighted mean transit depth, while dotted 
lines give l-cr uncertainty on this average value based on the scatter 
in the data. The numbers above each point indicate the orbital 
ephemeris. 



6.7. Eclipse depth variation 

Figure [U] shows the measured eclipse depth in units 
of the stellar flux at mid-eclipse. The weighted mean 
eclipse depth is 0.344 ± 0.0036% and the x^ fit to the 
eclipse depths with this mean value is 14.6 for 6 degrees of 
freedom. The errors on each eclipse depth vary between 
0.004-0.009%, while the scatter in the depths is 0.01%, 
so there is no detection of significant eclipse depth vari- 
ability. This scatter corresponds to 2.7% variation of the 
mid-eclipse planet brightness; this can be taken as an 
upper limit on the planet variability at 68% confidence. 



Although our eclipse depths do not exhibit signifi- 
cant variability, they only probe variability on timescales 
shorter than the baseline of the observations (roughly 
two years). To verify that our data do not show any 
time structure, we plot in Figure [T5] the change in eclipse 
depth against the time between observations, for each 
pair of observations {N x {N - l)/2 = 21 for N = 7). 
We add uncertainties in quadrature to estimate the un- 
certainty on the flux differences. The resulting locus is 
flat, showing that measured eclipse depth is uncorrelated 
with the time of the observation. Note that this sets a 
limit not only astrophysical variability scenarios, but also 
systematic errors. 
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Fig. 14. — Eclipse depths measured for seven planet eclipses. 
Horizontal solid line measures the average transit depth, while dot- 
ted lines give l-cr uncertainty on this average value. The numbers 
above each point indicate the orbital ephemeris. 
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Fig. 15. — For each pair of eclipse observations, we show here the 
change in eclipse depth as a percent, versus the time between the 
observations. If the eclipse depths showed time-correlation — due 
to either astrophysical variability on some characteristic timescale 
or detector systematics — this plot would show a rise. The flat 
distribution is consistent with Gaussian variations at the level of a 
few percent. 

This limit on the variation in eclipse depths is suffi- 
cient to rule ou t the predicted v ariati on computed for 
HD 189733b by iRauscher efaLl (|2008D . The mostex- 
treme prediction they make is for their rj — 0.05, U = 
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Fig. 16. — An exclusion plot for the covering fraction and temper- 
ature contrast of a putative storm on the day-side of HD 189733b. 
The top right corner of the plot is excluded at Scr. For comparison, 
Jupiter's Great Red Spot has a filling fraction of roughly 0.03, and 
a temperature 12 K cooler than the rest of the planet (asterisk). 



800 m s^^ model which has a standard deviation of ^ 8% 
in the day-side brightness, with the largest excursions of 
20%. The largest difference in brightness we see is 8%, 
while the scatter in planet brightness is 2.7%, so by both 
measures the observed variation is a factor of '^ 3 smaller 
than the predictions of this particular model. 

Other models pr edict smaller va riations in the day-side 
brightness, such as lShowman et al. (2009a ) who compute 
the 8 micron brightness variation for HD 189733b should 
be less than 1%. Our upper limits are consistent with this 
model, but unfortunately do not constrain the model due 
to our uncertainties that are larger than the predicted 
variations. 

With upper limits on both day and night-side vari- 
ability, it is worth asking which of these puts stronger 
constraints on the pla net's physical prop e rties. C onsider 
the simple model of iCowan and Agoll (|2010^ . which 
parametrizes the planet's day and night-side brightnesses 
in terms of the planet's Bond albedo and recirculation 
efficiency. One can treat brightness variability — of the 
day or night — as being due to changes in albedo and/or 
changes in recirculation efficiency, and compute how 
these affect the day and night side brightness. We find 
that the day side variability upper limit provides a better 
constraint on variation in the Bond albedo or recirula- 
tion efficiency than does the night side variability limit, 
which is ~ 6x larger than the day-side limit. 

Another possible origin of planet day-side variability 
are transient local variations in the surface brightness, 
for example due to large-scale "storms." We use a toy 
model where the planet's day side has a uniform tem- 
perature, T(j, except for a storm with covering fraction 
< / < 1 (the y-axis) and temperature difference AT 
from mean day-side temperature (the x-axis). Figure [T51 
shows exclusion limits on the largest putative storm that 
could form or dissipate without appearing in our data. 
The increasingly dark shades of gray denote areas of pa- 
rameter space e x cluded at la through 5a. According to 
iShowman et al.l ()2009bO . the radius of deformation for 
HD 189733b is 0.3 of the planetary radius, an order of 
magnitude larger than for Jupiter. The covering fraction 
for a typical storm on such a planet would be / = 0.1, for 



which we cannot rule out storms differing by 324 K 
confidence) from the average day-side temperature. For 
comparison, Jupiter's Great Red Spot has a filling frac- 
tion of roughly 0.03, and a temperature 12 K cooler than 
the rest of the planet. The bottom line is that our data 
rule out only the most extreme weather fluctuations on 
HD 189733b. 

6.8. System parameters 

Due to the high precision of our data and weak limb 
darkening in the infrared, we can considerably improve 
the determination of certain stellar and planet parame- 
ters for this system from our data. Since both the small 
inferred value of e cos uj and theoretical predictions in- 
dicate that e should be close to zero, we set e = in 
deriving the system parameters. The uncertainties on 
the stellar and planet parameters are computed for each 
transit or eclipse by computing the system parameters 
from the model p arameters fr om each simulation (us- 
ing the relations in lWinnll2010[ ). computing the standard 
deviation of the results from the simulations as an esti- 
mate of the errors on each parameter, and then taking a 
weighted mean of all transits/eclipses to obtain the final 
mean value of the best-fit parameters. 



TABLE 4 
Best fit system parameters 



Parameter 


Best fit 


units 


a/R, 


8.863±0.020 




fe/iJ. 


0.6631 ± 0.0023 




i 


85.710 ± 0.024 


deg 


ecosi..; 


0.000050± 0.000094 




Ul 


0.118 ± 0.010 




p* 


2.670 ± 0.017 


g cm-3 


Rp/R* 


0.155313 ± 0.000188 




Fp/F* 


0.3440 ± 0.0036 


% 


9p 


2145.9 ± 13.5 


cm s~'^ 


Pv 


0.943 ± 0.024 


gcm-3 



Table |4] presents the system parameters determined 
from all 14 transits and eclipses. We have focused 
on parameters that are most directly constrained from 
the photometry, which are either dimensionless, or have 
units of den s ity. Comp ared to the value s derived in 
[Torres et all (|2008[ ) and iPont et al.l (|2007D . the uncer- 
tainties on our values are smaller by a factor of 2-10. 
For the planet surface gravity, g^, we use the veloc- 
ity seni^ani£litude_^ = 200.56 ± 0.88 m s~^, derived 
bv iBoisse et al.l (|2009l ). and for the planet density, pp, 
w e use the ste l lar m ass M* = 0.806 ± O.O48M0 given 
in [Torres et al.l ()2008[ ). propagating the uncertainties as- 
suming they are uncorrelated and Gaussian. 

7. DISCUSSION AND CONCLUSIONS 

The analysis of fourteen transits and eclipses in this 
paper has made several improvements to the data reduc- 
tion and modeling; in particular, we have found a better 
function for fitting the detector ramp of IRAC Channel 
4, a double-exponential. The scatter in the residuals is 
approaching that of photon counting errors, similar to 
the precision achiev ed in other IRAC observations (e.g. 
iTodorov et al.ll20Tol) . but for a brighter source star, and 
the residuals show very little evidence of red noise. These 



16 



Agol et al. 



technical developments have allowed us to make a better 
correction for stellar variability, and have given us better 
constraints on the parameters of this system. 

As HD 189733 is the first planet system in which the 
phase variation has been measured at high significance, 
it provides some of the tightest constraints on the atmo- 
spheric dynamics of an extrasolar planet. Our observa- 
tions of an additional six transits and eclipses presented 
here allows us to place additional constraints on the lon- 
gitudinal brightness distribution of the planet at 8 mi- 
crons. In particular, we have improved the measurement 
of the correlation between the optical, {b + y)/2 band, 
and 8 micro n variations in t he sta r over the correlation 
measured in iKnutson et al.l ()2009l ). giving a correlation 
that agrees better with predictions of star spot models. 
This measured correlation allows us to derive a better 
correcti on for the steh ar variation during the observa- 
tion of IKnutson et al.l ([2007) , giving us a better mea- 
surement of the planet's phase function. In particular, 
we find that the peak planet flux at 8 microns occurs 
3.5 hours before secondar y eclipse, which is 1.2 hours be- 
fore the value derived in IKnutson et al.] ()2007D without 
correction for stellar variation; this is consistent within 
the errors given in that paper, which were dominated by 
the ramp correction. This measured phase function pre- 
dicts a 33 second delay of the secondary eclipse when fit 
with a uniform planet model, which is consistent with the 
38± 11 second delay that has been measured after cor- 
recting for light travel time across the system. This con- 
firms that the phase variation is indeed due to the planet, 
and gives a crude eclipse mapping of the planet detected 
the 3.5cr level, as first pointed out by IWilliams et al.l 
([2006). It is significant that — for the s ame high quality 
pho tometry — phase function mapping (i Cowan and Agol 
|200 8|) is more effective at locatin g the planet's primar y 
hot spot than eclipse mapping (jWilliams et al.l 120061 ) . 
This is because the duration of eclipse ingress or egress 
is shorter by a factor ~ Rp/{2Tra) than the planet's or- 
bital period, while the changes in brightness used by 
both techniques are comparable. The superior leverage of 
phase function mapping will become even more marked 
as interest shifts towards smaller planets in longer orbits. 
That said, the two mapping techniques suffer from dif- 
ferent degeneracies and different impacts of systematic 
errors and stellar variability, so when possible one will 
want to use both. We also confirm the phase variation 
by measuring the difference between the fiuxes at tran- 
sit and eclipse, and we find the night side is fainter by 
1.2±0.2 mmag, or about 64% of the brightness of the day 
side. All of these constraints are consistent with a model 
in which the gas circulating the planet has a radiative 
cooling timescale which is comparable to the advection 
timescale; we find Ti.ad/'''adv ^ 0.74 by fitting the phase 
function. 

The larger offset in the time of peak planet flux 
is also in better agre ement with the predictions of 
iShowman et al.l ()2009aD who foun d that to obtain agree- 
ment with the smaller offset of IKnutson et al.l (|2007D 
they required an inner boundary of their atmosphere 
that was rotating more slowly than synchronous rota- 
tion; instead of a sub-synchronous core, this may be in- 
dicative of slower wind spe eds due to magneti c drag near 
the 8 micron photosphere (jPerna et al.ll2010l) . The sub- 
synchronously rotating and 5x solar abundance models 



of lShowman et al.l ()2009aD predict a peak brightness at 8 
microns which is 20-30 degrees before secondary eclipse, 
which agrees well with our new estimate of 27 degrees. 
The same models also predicts a day-side brightness 
(mid-transit) which is 0.33-0.35% of the star's brightness, 
consistent with our measured value of 0.344±0.004%. 
The night side brightness at 8 micron predicted by the 
models is 0.17-0.18% of the stellar brightness, which is 
consistent with our measured value of 0.22±0.05%. Their 
models also predict very small variations in the secondary 
eclipse depth of less than 1%, which is consistent with our 
upper limit of 2.7%. The lack of variation of the atmo- 
sphere indicates that the assumptions used in creating 
longitudinal maps of planets from phase functions are 
likely vahd ()Cowan anid Agolll2008( ). 

The time delay for the secondary eclipse can be com- 
pletely accounted for by the light travel time of the sys- 
tem and delay of ingress and egress due to a hotspot on 
the planet which is offset longitudinally. Consequently 
there is no evidence for a non-zero ecosw, and we can 
place a limit of ecosw = 0.00005 ± 0.00009. If the orbit 
of this planet is nearly circular, which the small value 
of ecosoj would indicate, then the interior is likely also 
synchronously rotating , which seems to agree with the 
IShowman et al.l ()2009af) predictions for the phase func- 
tion. 

We have detected 8 /xm limb-darkening of the star at 
high significance, ~ IOct, which agrees with predictions 
of stellar atmosphere models. However, the individual 
transits vary in depth, which we hypothesize may be due 
to variation of the stellar surface brightness that is oc- 
culted by the planet. This is not surprising given the 
strong optical variations of this star which indicate a sig- 
nificant presence of star spots. This variation needs to 
be accounted for in creating spectral absorption profiles 
of transiting planets. If the data taken are non- simul- 
taneous, the variation in stellar surface brightness could 
affect the inferred depth of transit differently at different 
wavelengths, leading to systematic errors in comparison 
to model predictions; even simultaneous data might be 
affected by the star spot color. This is a stronger ef- 
fect at shorter wavelengths; for example, the contrast 
in surface brightness of 4000 K star spots in the IRAC 
Channel 1 (3.6 /im) and Channel 2 (4.5 /im) should be 
20-40% higher than for Channel 4 for this star; thus fiuc- 
tuations in transit depth could approach 2% in these 
bands. In addition, this limits the possibility of con- 
straining t he variations in transit depth due to planet 
oblateness ([Carter and Winnll2010[ ). 

Due to our highly precise transit times spaced over 
a wide range in time, the ephemeris we derive is one 
of the most precise for any transiting planet. The 
high precision is due to the weak limb-darkening, sta- 
ble instrument (thanks to the Earth-trailing orbit of 
Spitzer which leads to stable thermal properties and 
no occultation of targets by the Earth, as occurs with 
the Hubble Space Telescope), allowing a 3-second pre- 
cision for transit times. Our ephemeris has a precision 
that is > 10 times better for the period than that re- 
ported in in IPont et all ([20071 ). and agrees with their 
reported period within ~2.8i7: their period is longer by 
0.46±0.17 sec. Our ephemeris predicts times of transit 
that are -3.3±5.0, 3.5±5.0, and 12.6zb3 .5 seconds after 
their three transit times in ,Pont et al.l ([2007, ). We de- 
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tect no strong evidence for transit timing variations in 
our data, and we estimated from analytic formulae the 
upper mass limits on the presence of companion plan- 
ets in this system , impr oving upon the limits placed by 
iMiller-RiccieTall ([200^ and iHrudkova etHI (poT o'l by 
a factor of ^ 10. Theories of the evolution of short- 
period planets due to tidal effects and interaction with 
turbulence in the protoplanetary disk indicate that they 
should evolve out of mean-motion resonance, so the lack 
of detected transit-timing variations may not be surpris- 
ing, especially for interior p erturbing planets (Fab rycky 
2009; Adams' ct al. 2008; Te rmiem and Papaloizoull2007l: 
Papaloizou and Teraucm..2010l ). 



In sum, the excellent stability of the Spitzer Space Tele- 
scope, and in particular Channel 4 of the IRAC camera, 
has enabled near photon-limited photometric errors, and 
sub-mmag variations over a period 1.6 years. This has 
enabled a better calibration of the contribution of stellar 
variability at 8 micron, which has allowed us to mea- 



sure the night-side planet brightness, and has shifted our 
estimate of the peak of the phase function. 
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